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SUMMARY 

A very general method for calculating compressible three-dimensional laminar and 
turbulent boundary layers on arbitrary wings is described. The method utilizes a non- 
orthogbnal coordinate system for the boundary-layer calculations and Includes a geometry 
package that represents the wing analytically. In the calculations all the geometric 
parameters of the coordinate system are accounted for. The Re)molds shear -stress 
terms are modeled by an eddy-viscosity formulation developed by Cebeci. The govern- 
ing equations are solved by a very efficient two-point finite -difference method used ear- 
lier by Keller and Cebeci for two-dimensional flows and later by Cebeci for three- 
dimensional flows. 

Preliminary results for a swept wing look very encouraging. A typical computation 
time (CPU) for one surface of the wing which roughly consists of 30 z -stations and 20 
stations with 30 ? 7 -points across the boimdary layer is a little over 30 sec on an 
IBM 370/165 computer. 

INTRODUCTION 

The development of an efficient and accurate method to compute three-dimensional 
boundary layers on wings of arbitrary shape requires: 

(1) The velocity distribution at the boundary -layer edge 

(2) A convenient coordinate system 

(3) A model for the Reynolds stresses 

(4) A numerical method to solve the governing equations 

The velocity distribution must be obtained from the pressure distribution. In gen- 
eral, the pressure distribution can be obtained either theoretically or experimentally. 
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When obtained theoretically, the velocity components in the streamwise and spanwise 
directions can be' calculated without too much difficulty and thus satisfy the first require- 
ment. When the pressure distribution is obtained experimentally, the calculation of the ■ 
velocity components ’is rather difficult; Certain approximations must be made to get the 
velocity distribution from the experimental pressure distribution. In the section "Govr . 
erning Equations,!' the difficulties and the procedure used to calculate the velocity com- 
ponents from the experimental pressure distribution are discussed. 

‘ In selecting a coordinate system for Uie boundary -l^ayer calculations, aji important. __ 
point to consider is thatjthe coordinate system should be calculated only once for each 
geometrical configuration. This rules out the streamline coordinate system since for . 
each angle of attack the streamlines must be calculated repeatedly. Another important 
point to consider is dictated by utility. The measured or calculated external velocity 
distributions are usually given in planes containing the local chord line. Hence it is 
natural to select one surface coordinate in planes parallel to the defining sections. The 
other surface coordinate may be lines either orthogonal or nonorthpgonal to that coordi- 
nate line. However, the selection of an orthogonal system causes a number of inconve- 
niences together with lengthy interpolation procedures. As a result, a nonorthogonal 
coordinate system appears to be the most convenient system with which to perform the 
boundary -layer calculations as discussed in detail in. the section "Coordinate System." 

For turbulent flows the governing boundary-layer equations contain the Reynolds 
stress terms which require closure assumptions such as mixing-length, eddy-viscosity 
concepts or "higher order turbulence" models. Although the latter have the potential to 
compute more complicated turbulent flows, mixing -length, eddy-viscosity approaches 
have proven to yield quite satisfactory results for boundary-layer flows. (See refs. 1 to 4.) 
The use of higher order turbulence models also increases the complexity of already com- 
plex equations leading to high computation times. Furthermore, for compressible flows 
their accuracy may not be as good as the simple mixing -let^h, eddy -viscosity methods. 
For this reason, in our study the Re 3 molds stresses are modeled by using ai^ accurate ; 
eddy -viscosity formulation (see the section "Turbulence Model") developed by Cebeci. 

(See refs. 3 and 4.) 

. When physical coordinates are used, the solutions of the governing boundary^layer 
equations are quite sensitive to the spacings in the streamwise direction (x) and to the 
spanwise direction (z) and require a large number of x- and z -stations. In calculations 
such as the ones considered here where the computation time and storage become impor- 
tant, it is necessary to remove the sensitivity to Ax- and Az -spacings. This can be done 
by expressii^ and solving the governing equations in transformed coordinates. There-, 
fore, in the section "Transformation of the Governing Equations," a convenient transfor- 
mation to ejq)ress the boundary -layer equations in terms of transformed variables is 
considered. 
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In the section "Numerical Method," the solution of the governing equations the 
Cebecl -Keller Box method is discussed. This is a very efficient two-point finite- 
difference method develcq}ed by H. B. Keller and applied to the boundary-layer equations . 
by Keller and Cebeci. (See, for example,, refs. 5 and 6.) 

In the section "Results and Discussion," results for a planar turbulent boundary- 
layer flow approaching a three-dimensional obstacle and results for a swept wing are 
presented. Finally in the section "Futiire Work," additional work that needs to be done 
in order to develop a complete design tool for computing, the flow field past an arbitrary 
wing is discussed. 

SYMBOLS 

A Van Driest damping length, 26(i'/u7)(p/p^)^/^ 

b =c(i+€+,) ; , . . : : 

C = P4/(PePe) . \ 

Cp pressure coefficient, 2(p -p„)/(Pi,Uoo2) 

c the ratio Pq/p; local chord 

Cf skin -friction coefficient 

E total enthalpy ratio, H/Hg 

f transformed vector potential for ^ 

G = E' 

g transformed vector potential for ^ 

g' = W/Wg 

H total enthalpy 
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h,hj,h2 metric coefficients 

T,J,? imlt vectors In x-, y-, and z-dlrectlons of Cartesian coordinate system in 
which wing Is defined 

2 curvature vec^r of coordinate line 

/ geodesic curvatures 

K]^2>^21 geometric parameters 

L modified mixing length ' ' ’ 

Moo free-stream Mach number 

Npr Prandtl number 

n unit vector normal to the surface ' 

P parameter denoting either coordinate, <t> or y; point 

• M^IO parameters In transformed differential equations 

p static pressure 

p^ total pressure 

= vl^8i/vq, Reynolds number 

r see figure 7 

r position vector for point on surface, ^,y,z) 

S stagnation point 

8 distance along a streamline 

/ 
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Si 

t 


curve length along x -coordinate line 
’ time ■ 



t unit tangent vector along coordinate line . 

u,v,w,t dependent variables in first-order, transformed differential equations^ f, 
r, g’, and g” 

u see figure 7 

Ug total velocity (ut evaluated at the edge) 


ut total or resultant velocity 

Ur friction velocity, 

-A- '' 

V velocity normal to surface in physical differential equations , 

x,z,y independent coordinates in boundary-layer equations 
x,y,z Cartesian coordinate system used for wing definition 

a local geometric angle of attack of wing section chord lines 

p = Sin-l(Up/Ug) _ ^ , fn: 


y ratio of specific heats, y = 1.4 ' 

e(or eddy viscosity Md eddy conductivity, respectively '• 

7} transformed coordinate normal to surface 

angle in tangent plane between x- and z -coordinate lilies 

\ local sweep angle, measured between plane normal to free -stream velocity 

vector and z -coordinate line 

t ' ■ ^ 

j 
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molecular viscosity 


M 

parameters in transformed energy equation 

V kinematic viscosity 

p density 

%w resultant wall shear stress 

shear stresses 

0 stretching variable defined in figure 3 

4,0 two-component vector potentials, equation (56) 

Subscripts: 

e outer edge 

g geodesic 

1 input stations 

y,n . indices . - 

In, out inner and outer regions for eddy viscosity 

le leading edge 

p pivotal points, l.e., stations at which boun^ry layer is computed 

t wing tip 

te trailing edge 

w wall 

\ 
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free -stream conditions 

c 

Primes denote differentiation with respect to 77. 

GOVERNING EQUATIONS 


The governing boundary -layer equations for a non orthogonal coordinate system are. 
given in references 7 and 8. With a slight change of notation for compressible laminar 
and turbvilent flows, they are given by - , 

Continuity equation: 

-|^(puh2 sin 0) + ^(pwhj sin 0) + ^(^hih2 sin 0) = 0 (1) 

X -momentum equation: ' ■ ; ■ • ' ^ < i 

p ^ + p ^ + pv — - pKiu^ cot 0 + pKow? CSC 0 + pKiouw . , . 

hi ax h2 az ay ^ 1 


ax 


= _ csc.^l jk ■ cot 0 CSC 0 ap . 3 /„ ^ 
h' 

z -momentum equation: 


h2 


dz ay^ ay ) 


cot 0 + pKiu2 csc 0 + PK21UW 


_ cot 0 CSC 0 9p csc^ 0 ap a 


9 CSC 0 op csc^ 0 £P + ^ /m aw 
hi ax ‘ h 2 az ay\^ ay " j 


Energy equation: 


hi ax h2 az 


a 


— + u /l 

1 ' 



ay 

Npr 

8y *‘^1, 

'Npr; 

)ay' 

\2j 




(4) 


where pv = py + p'v' and hi. and h2 are metric coefficients. The latter ^are func- 
tions of X and z, that is, 

hi = hi(x,z) h2 = h2(x,z) (5)'- 


Also, 0 represents the angle between the coordinate lines x and z. F<y ^ orthog- . 
onal system 0 = ff/2. The parameters Ki and K2 are known as the geodesic curva- 
tures of the curves z = Constant and x = Constant, respectively. They are-given by 
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Ki i COS 6 ) - 

^ hjh 2 sin 0|_9x' 2 / 


— - -^(hi cos 0) - 
in ^ ' 


hxh 2 sin 


The parameters Kj2 and K 21 are defined by 


Ki2 = ^ 


sin 9 




''a 1 = 3 ih -^''2 * it ill * ® (’■'1 ^ » 


^2 ^z 


hi ax 


The total velocity within the boundary layer uj is given by 

uj = (u^ + + 2uw cos 9 ) ^ (8) 

One obvious procedure to calculate the velocity components Ug and We from the 
given pressure distribution is to evaluate equations (2) and (3) at the edge of the boundary 
layer. This gives 

£iie ^ ^ _ K ug2 cot 0 + K2We2 esc 9 + Ki 2 UeWe = - ^ + cot0_c^ ^ 

hi ax h2 az 1 e 2 6 12 e e 


Ue awe . We awe 
h^’aF'^ 


K 2 We^ cot 0 + KjUe^ esc 0 + KoiUeWe = £ 9 . t Q esc 0 |P _ e sj;^ 9 |P 

" ■ ^ h-iD. aX h«n oZ 


Equations (9)' and (10), which may be expressed in the form 
He 9«e . We 9ue _ , ^ ^ 

hi IF * hi- IT ' 


= (12 

constitute a system of first-order quasi -linear partial differential equations in Ug and 
Wg. The differential relationships for these variables are 
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(13) 


dUg = 


9ue 

ax 


dx + 


aue 

dZ 


dz 


dwg = dx + dz . (14) 

® ax 3z 

If we let s denote the distance along a streamline, and Ug the total velocity (ut eval- 
uated at the edge), that is. 


Us = (u«2 + Wg2 + 2ugWg cos 0 ) 


1/2 


equations (13) and (14) can be expressed in the form 
Us 


dug _ Ug aug wg aug 


ds 


hj 9x 


h2 9z 


(15) 


(16) 


and 

u dWg _ Ug 9wg ^ Wg awg 
^ ds hj 9x h 2 9z 

by noting that 

Us=f = Ws = h2| 

Comparison of equations (11) and (12) with equations (16) and (17) gives 
= X dWg _ G 

ds Ug dg Ug 


(17) 


(18) 


(19) 


In addition, we have the following relationships 

^ - dz _ Wg _ 

ds hjUg ds h2Ug 

The system of four first-order differential equations (eqs. (19) and (20)) allows one to 
calculate the variation of ug, Wg, x, and z along a streamline. In principle, these 
eqviations can be solved as an initial-value problem. However, it can be shovm that the 
system of differential equations (11) and (12) has characteristics which are identical to the 
inviscid streamlines. As a result, the initial-value problem cannot be started from lines 
which are streamlines. Thus, with initial points on the stagnation line or in the plane of 
symmetry, the solution is quite difficult except for the initial lines themselves. To obtain 
the solutions over the entire surface, the initial values of Ug and Wg must be known 
along a line which is not a streamline itself. However, this information is not available 
in general. A satisfactory solution requires considerable study. In this study approxi- 
mate methods are used. The simple sweep theory is known to give reasonable answers 
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when applied to regions of high-aspect -ratio wings that are outside the influence of root '■ 
and tip effects. In the absence of spanwise pressure gradients, this approximation is , 
almost exact. Thus, for weak spanwise pressure gradients, we can obtain the velocity., ■ 
components on” the midportion of a swept wing with reasonable accuracy by using the ■ ; , . 
sweep theory. In regions of root and tip influence, the simple sweep theory with a cor- 
rection to the sweep angle is applied. The procedure is explained below. 

Consider the velocity vector in the tangent plane at a point P on the wing. (See 
fig.. 1.)- Thebasic assumption for the simple sweep theory is that the velocity ‘component 
xip in the z -direction is given by: 

Up = Uoo sin X. (21) 

The sweep angle X represents the angle between the spanwise direction and the z- 
coordinate line through the point P. The parallelogram addition of vector coniponents 
yields 

cosf . 

Uoo U«, sin 0 ^ ’ 


We _ Us sin sin 0 - cos 9 cos 
Uoo Uoo sin 0 

where sin ^ SL^ : . 

Elimination of /3 from equations (22) and (23) yields 
^ ^ \lius/u^)^ - sin2 

^ sin.0 _ .. . — - 

■^ = sin X - ^ cos 0 

Uoo Uoo 

The total velocity ratio ug/uoo is calculated from 



1 - 


m- 


= 1 + 


Pj,2\ 2 


; (y-l)/y 


CpM«.‘ 


(y - 1 )Moo2 


(23) 


(24) 

(25) 

(26) 


with Pj j and p^ g denoting the va,lues of total pressure before and after the shock, 
respectively, and ’Cp is the pressure coefficient, Cp = (p - Poo)/(l/2pUoo2). Equa- 
tion (26) is valid for an adiabatic flow through a shock wave, but since the total pressure 
ratio across the shock is seldom known, its effect will be neglected. This approximation 


50 



introduces only an error of a few percent into the velocity calculations because the total 
pressure jump across a swept shock is small even for free -stream Mach numbers 
approaching unity. The total pressure ratio must also remain close to one for the first- 
order boimdary-layer theory to be valid in'front of arid behind the shock wave. . ' 


Equations (24) to (26) are approximately valid for the root and tip regions if the 
local sweep angle X is replaced by an effective sweep angle Xgff 

^eff = ^ - ?r^i^ 

> (27) 

^eff = ^ - FtXt J 


where Xf and X^ denote the root and tip sweep angle for the given z -coordinate line 
and Ff and Fj are the spanwise interpolation factors for the root and tip, respectively. 
These parameters are shown schematically in reference 9 as a function of nondimensional 
spanwise distance in terms of root or tip chord. 


COORDINATE SYSTEM 


The wing is defined in the x,y,z coordinate system. Here, the x-axis is in the 
direction of the airplane's longitudinal axis and the y-axis is in the spanwise direction. 

It is assumed that the wing is defined by a number of airfoil sections in the planes 
y a Constant, which involve the specification of Z[ and x^ for constant values of yj. 

It is also assumed that the pivotal points along the chordwise direction (x/c)p where c 
denotes local chord are given, as are the spanwise stations yp where the boundary- 
layer calculations are to be made. These parameters are shown schematically iri 
figure 2. 


The defining airfoils are usually given by n pairs of values of xj and zj. But 
because all aerodynamic data related to airfoils are customarily given in terms of frac- 
tion of the total chord c, the input data are converted to an xz -coordinate system (see 
fig. 3) based on the local chord (maximum length line). The relationships between x, 

X, z, and z are 


I = - *le) cos o - (z - zje) sin oj 

S = - *le) Of + (z - zie) cos oj 


(28) 


(29) 


where 


c 




= (xte -xie)^ + (zte - ^le)^ 


(30) 
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( 31 ) 


a = -tan“l ^ 

^le-^e 

The subscripts le and te refer to the points at the leading and trailing edges. They 
must be specified. 

The curves yj = Constant and the curves connecting the points x/c = Constant 
on all the defining airfoils form a convenient coordinate system. However, the movement 
of the stagnation point S with.angle of attack gives, rise to- ambiguity. — ^For example; 
the same x/c value may correspond to two z/c values on a given section. To avoid 
this problem, another variable <t> defined by 

. = |(1 - cos 0) (32) 

is 'introduced. Here, 0 = 0 corresponds to the leading edge, 0 = tt corresponds to the 
trailing edge. The value of 0 is positive for the upper surface. Qi the lower surface 
0 is negative. 

Other Possible Coordinate Systems 

•’ In this section, other possible coordinate systems are discussed. Because of - 
impracticalitles with these systems, the nonorthogonal coordinate system is the most 
convenient to perform the boundary -layer calculations for wing surfaces. 

As pointed out in the Introduction, one surface coordinate must be chosen in planes 
parallel to the defining sections. Consider an orthogonal system in which the orthogonals 
are construrted between the intersections of planes parallel to the defining sections and 
the wing surface. Trial calculations showed that orthogonals started from the wing -root 
congregate. at the leading edge, leaving large portions of the wing uncovered. (See fi'g. 4.) 
This is especially true for a wing with a sharp trailing ec^e. A rounded trailing edge 
rectiiies'the situation somewhat but there is still a large area of the wing where the 
orthogonals are sparse. • . , : 

The orthc^onal coordinate system in figure 4 was constructed with the polar angle 
0 = X 2 at the root section as the other surface coordinate. As is seen from this figure, 
there are’ computational difficulties at the trailing edge. To show this, consider figure 5, 
in which the surface coordinates xj and X 2 are obtained by extending the surface cov- 
erage with the dashed lines. Here, AA" is the stagnation line, AB is the root section, 
and D is a point on the trailing edge. Starting from the initial lines, the boundary layer 
can be calculated along the line BC" including the root chord. However, the point D 
cannot be obtained in a straightforward manner. This is also true for the rest of the 
trailing -edge points D', 0', and D'". Because of the difficulty in' calculating these 
trailing-edge points, the orthc^onal coordinate system is not practical. 


52 



Another possible coordinate system can be obtained by representing the wing by one 
or more separate conical surfaces. Figure 6 shows such a representation. Here, the 
wii^ panels ABDC and CDFE form two conical surfaces with apexes at P and Q, 
respectively. The shape of the panel ABDC and the coordinate system in the developed 
plane are shown in figure 7. The initial lines are AC and AB. Line AC is the stag- 
nation line and AB is the wing -fuselage junction. Calculations can be started at corner 
A. A linear coordinate transformation can be used to avoid marching into the negative r 
direction. Such a coordinate system without taking the thickness into account (this 
amounts to representing wing sections by flat plates) was used by Nash and Scruggs 
(ref. 10). The disadvantage of this coordinate system is the difficulty of doing calcula- 
tions in the overlz^ region. 

' . • f'.. '■ 

Present Coordinate System 

The most convenient coordinate system on the wing surface and the one used in this 
study is a nonorthogonal coordinate system given by the lines y = Constant and 
<f> = Constant. The new independent variables <p and y are selected to correspond to-, 
the independent boundary -layer parameters x and z, respectively, in equations (1) to 
(4). Before the boundary-layer calculations are performed, it is necessary to decide on 
the surface locations for which the boundary -layer solutions will be output. The best 
method is a chordwise point distribution in terms of percent chord. The information can 
then be converted to give the <^>p -values for the pivotal points. As is likely to h^pen, the 
points on the wing defining sections will not correspond to the pivotal points for the ' 
boundary -layer calculations. Thus interpolation is necessary. At each spanwise defining 
station, the corresponding to the input data can be found by using equations (28) to 
(32). Next, Xj versus and Zj versus are curve fitted with cubic spline func- 
tions and Xpj and Zpj are interpolated for at each sp^wise station. Then the vXpj 
and are spline fitted versus yj for each (p^ and are interpolated for Xp and 

% at yp. 


Calculation of the Geometric Parameters of the Coordinate System 

Once the coordinate system is selected, it is necessary to calculate its geometric 
parameters, namely, the metric coefficients hj and h 2 and Kj and K 2 which 
appear in the governing boundary -layer equations. These are calculated by the procedure 
described below. 

The metric coefficient along one curve in space is given by 




with P denoting a parameter. For P = <f> along the curves y = Constant, equation (33) 
can be written for hj as 


h.2 = 


y 


(34) 


Similarly, for P = y along the curves 0 = Constant 

/.sTr\2 


h22 = i + f|ir +« 

A, 


'0 


(35) 


The derivatives in equations (34) and (35), namely (ax/30)y, (3z/a0)y, (ax/3y)^, and 
(9z/9y)0, can be obtained as byproducts of spline -fitting the points along the chordwise 
and spanwise directions at the pivotal points. 


The unit tangent vector t along a curve is given by 


T = ^ 


_ 1 dr 


ds dP ds/dP h dP 
The unit tangent vector t j along the curve y = Constant is 




where i, j,and k are unit vectors in the coordinate directions x, y, and z, 
respectively. Thd unit tangent vector T 2 along the curve 0 = Constant is 


/9x\ J ^ 


\90/y \d(f. 

7y J 


(36) 


(37) 


to = ^ 
^ ^2 


T + T + k 

\9y/0 \9y/0 


(38) 


The angle between the coordinate lines is then 


fdxWdx\ , /azVszN 

T r \90A9y/ \90A9y/ 

= ti . t2 = — 


(39) 


cos 6 = . 

hih2 

The curvature of a curve in space is given by 

g - dT _ dt 1 _ 1 dt 


(40) 


ds dP ds/dP h dP 

The geodesic or tangential curvature Kg of a curve on the surface can be obtained from 
Kg=(T><n)-K (41) 

i . . 

Here, n is the vector normal to the surface which by definition is 

if sin 0 - 1 1 X T 2 (42) 
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or 


n = 


hj^h 2 sin 9 


- 

/dx 8Z 

8x 9 z\T , /a 

i-'lk 

V90/ 

\90 8y 

9y 90y^ \a 

0/ 


(43) 


With the use of equation (36), equation (40) can be written as 

K = A ^ X dr dh ^ x/d% T ^ d% ^ 1 (dS r , ^ T . dz irVd%c cS 


h2 (jp2 ^3 dP dP h2\dp2 ^ ^ 


dP' 


] + 


k -X 


dP^ 


4\dP 


i + 


dP 


dP 


J t W U A UA 


dP 


2 dP 


+ ^ 


dp2dp - dp2dP/ 

The geodesic curvature for a curve y = Constant is 

Kg^ = -(t 1 X n) • Ki (45) 

The minus sign on the right-hand side of equation (45) is introduced to obtain 

Kg^ = - (l/hih 2 ) (3hj/az) in the case of an orthogonal coordinate system. With as 

the parameter, the expression for Kj is 


T? 1 /d^X rA 

hi^\30^. 


Ill 

90' 


k - 


1 


/9X ![■ 9^ g^\ /9%C 8x ^ 9% 

\80 80 /l802 90 302 


3zj 
302 901 


(46) 


Substituting equations (37), (43), and (46) into equation (45) gives, after simplifications 

l’ (dx dz 9X dz\ fd^ 9z 9% 9xV 


hi'^h2 sine\^^ ®y 9<^/\902 90 8.02 90y 

The geodesic curvature Kg^ for a curve 0 = Constant is given by 

•%2 = (‘2 X if) • K2 

With y as the parameter, the expression for K£ is 

h2^\9y^ 8 v 2 / ho^\®y ®y Aav^ ^y 8v2 ®y/ 


(47) 


(48) 


(49) 


h2^ 


The expression for the geodesic curvature Kg^ is obtained by substitution of equa- 
tions (38), (43), and (49) into equation (48): 


%2 = 


hjh 2 ^ sin 9 


/9X 8Z 9X 9z\/ 9^ 9Z 9f|z 9z 

1,80 8y ■ 8y 80yL^2 df ' gy2 8y/ ' .lg 3^2 80 ' gy2 80 


^8% 82L + 8^z 8z^ 


(50) 
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The second -order partial derivatives appearing in equations (47) and (50) are also 
obtained as a byproduct of the spline -fitting technique. In terms of parameters appear- 
ii^ in boundary-layer equations, we set Kj = and K 2 = Kg^. In addition to 6 , 
hl> ^2» ^1> bovmdary -layer equations contain K 12 and K 21 which are 

functions of the previously mentioned parameters. Also, the partial derivatives, 30/90 
and 30/3y, are contained in the boundary -layer equations and obtained by spline-fittii^ 

0 versus 4 > and y. 


TURBULENCE MODEL 


The solution of the system of equations (1) to (4) requires closure assumptions for 
the Reynolds stresses, ’ -puV, -pwV, and -pv'H'. This can be done by a number of 
approaches. One approach is to use simple eddy -viscosity and mixing -length formulas 
for the Reynolds stresses. This method, also called the mean-field method, has been 
used by Cebeci and Smith (ref. 11), Bushnell and Beckwith (ref. 1), and Harris (ref. 2) as 
well as several others. Another approach is to use expressions that consider the rate of 
change of the Reynolds stresses in the governing equations. This method, called 
transport -equation method, has been used by Bradshaw (ref. 12), Donaldson and SuUivan 
(ref. 13), Hanjalic and Launder (ref. 14), and several others. In reference 15, Bradshaw 
presents an excellent discussion of both these methods. 

For low -speed flows, both approaches work equally well. For high-speed flows, 
however, the mean-field method seems to be slightly better than the transport -equation 
method, chiefly because of the inadequate closure assumption accounting for the mean 
compression or dilatation effect. However, a recent report by Bradshaw (ref. 16) seems 
to improve substantially the predictions of his method for compressible flows. In either 
case, equations (1) to (4) are already quite difficult to solve, and there is no need to 
increase the computation time by using higher order turbulence models. For this reason, 
an eddy -viscosity formulation developed by Cebeci (refs. 3 and 4) is used in this study. 
According to this formulation, the boundary layer is divided into two regions, called inner 
and outer regions, and eddy-viscosity formulas are defined separately in each region. 

For a nonorthogonal system (assuming no mass transfer), the inner eddy viscosity 
is defined by 



where 

L = 0.4y[l - exp (-y/A)] 


(51) 


(52) 
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The total shear stress evaluated at the wall is 

v 2 / q ™\2 


"^t.w = 


'Mr +f|wr + 2 cosef'M w 
Ww ry^w ry/wry/wi 


1/2 


The outer eddy viscosity is defined by the formula 


^m,out “ 


ppo 

\ (“t,e - “t)<iy| 


where 


'^t,e = (“e^ + + ZugWe cos e) 

uj = (u^ + w^ + 2uw cos 0 ) ^ 
and a = 0.0168. 

TRANSFORMATION OF THE GOVERNING EQUATIONS 
Boundary-Layer Equations 

Two-component vector potentials \f/ and ^ are defined such that 
puh 2 sin 0 = 1^ 


pwhj sin ® 
^hjh 2 sin 0 = -I 


The following transformations are also defined 

X = X 

z = z 

A I ^e A 

= TTT .^cJ p dy 


PePe®l/ 

= (PeMeUeSi)^/^h 2 f(x,z,?;) sin 0 

$ = (PePe'ieSi)^/^ ^ hj g{x,z,T]) sin 0 


(53) 

(54) 

(55a) 

(55b) 


(56) 


(57a) 

(57b) 

(57c) 

(57d), 

(57e) 
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where 




dx 


(58) 


Substituting equations (56) and (57) into equations (2) to (4), after considerable algebra, 
gives the following: 

X -momentum equation: 


(bf’)' +“Piff’ + P2[c -'(f’)2j+ Ps(c - f'g’) + Pgf’g + Pa[c - (g’)2] = xPio[f’ f 1^ 


(59) 


z -momentum equation: 

(bg")’ + Pifg" + P4(c f’g’) + P3[c - (g’)2] + Pggg” + Pg[c - (f’)2] = xPjo 


f’ _ e" M. 
ax ^ ax 


+ P. 


.fg* . gH IS) 
az ^ dzj 


Energy equation: 

(|i lE’)’ + M2E’ + = xPio 


tt aE pt af , p_/^« aE pr 9g\ 
f ^-E ^ + P7(g a^--® W 


(6b) 

(61) 


Here, primes denote differentiation with respect to t ] and 

P, =1 + ^ l_(p u ) - Si /ki cot 9 ^12 cos ^ + ^21 

2 2uehj ax 2pefiehi ax^®^©^ 1 gj^^2 q 


P2- 
P3 = 



®1 

2 2ueh]^ 

Si 

aue 

uehi 

ax 

Si 

awe 

Ueh2 

3z 

®1 

8We 


ax 


- Kjsj cot 0 


Up 


K21S1 


(62a) 

(62b) 

(62c) 

(62d) 
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Ps 




Si 3vie\ 

Ueh2 J 


p - We ^1 9We _ Si 3ue Sj j_ v JL £fl 


[weh2 2ueh2 2pg.fi e^2 


2h2 


dz 


- si(K 2 cot 9 - 


K\2 + ^21 ^®s o' 

sin2 e 


P7=IilZ[e 

7 h2Ue 


Pa = (sj) K2S1 CSC 9 


(62e) 


(62f) 

(62g) 

(62h) 


Pg = ^ Kisi CSC 9 


(621) 


and 



(62j) 


Npr' 

\ 1 

Npr,t 

/Npr 


M2 = Pi« + Peg 
2 


M3 




vr + ^ g'g" + cos 9 ^(g’f" + f'g") 
ue^ “e 




c -_£H_ 

Pef^e 


b = C(1 + eJn) 


c = ?e 


(63a) ' 

(63b) 

(63c) 

(63d) 

(63e) 

(63f) 

(63g) 
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In the preceding equations, eddy -viscosity and eddy-conductivity concepts have been used 
in order to satisfy the closure conditions for the Reynolds stresses. .They are defined by 

-puV = pem 1^ -pV^ = pem -p^ = P^h (64) 

The turbiilent Prandtl number Npj. j and the dimensionless transport coefficients are 
defined by 


-Npr,l-=# - 






(65) 


Equations (59) to (61) are subject to the following boundary conditions: 

t] = 0 f = g = 0 f = g' = 0 E' ~ 0 (adiabatic wall)! 


V ^ rjc 


V ~ 1 


r’ -* 1 


E - 1 


( 66 ) 


Eddy-Viscosity Equations 

The eddy -viscosity formulas given by equations (51) to (55) can also be transformed 


and e?q)ressed as 


'm,m ~ ~pr 




V2 


(67) 




1 + —(2 COS 6 + -^ 


Ue 


.1 1/2 


u, 




1/^ 




}dTj 

/ ) 


where 




C JQP 


(Iw)^ + § g'w(2fw cos 0 + g g'^^ 


1/4 


(68) 


(69) 


NUMERICAL METHOD 

The Cebeci -Keller Box method is used to solve the governing boundary -layer equa- 
tions given by equations (59) to (61), This is a two-point finite -difference method devel- 
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(^ed by H. B. Keller (ref. 17) and applied to the boundary -layer equations by Keller and 
Cebeci (re’fs. 5 and 6). The method is discussed in detail in references 6 and 11. For 
this reason only a brief description will be given here. 

C^e of the basic ideas of this method is to write the governing system of equations 
in the form of a first-order system. Thus, derivatives of some quantities with respect 
to the "normal" variable must be introduced as new unknown fimctions. Derivatives with 
respect to all other variables occur only to first order as a consequence of the boundary - 
layer approximations. With the resulting first-order system and an arbitrary rectangular 
net, centered difference quotients and averages at the midpoints.of net rectangles and net 
segments are used, as required, to get 0(h2) accurate finite -difference equations. 

This method is unconditionally stable; however, the equations are highly implicit 
and nonlinear. Newton's method is employed to solve them. In order to do this with an 
efficient and stable computational scheme, a block -tridiagonal factorization scheme is 
used. 

Numerical Formulation of the Momentum Equations 

New dependent variables u(x,z,?7), v(x,z,?7), w(x,z,t;), and t{x,z,rj) are intro- 
duced, so that equations (59) and (60) can be written as 

(bv)' + P^fv + P2(c - u^) + P5(c - uw) + P0gv + Pg(c - w^) = xPjq 


u 


8x 


9l 

8X 


+ P-; 


- V ^ 

V 8z dzj 


(70a) 


(bt)' + Pjft + P4(c - uw) + P3(c - w^) + P0gt + Pg(c - u^) = xPjg 


u|2i.t|l 
8x 8x 


+ P, 


'\ 8Z 8Z> 


(70b) 


f =u 
u' = V 
g' = w 
w' = t 


For the net cube shown in figure 8, the net points are 

Xo = 0 Xn = Xn_i + kn (n = 1, 2, , . ., N) 


(70c) 

(70d) 

(70e) 

(70f) 

(71a) 
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Zo = 0 Zi = Zi_i + Tj , (i = 1, 2, . . I) (71b) 

T?o = 0 7Jj=7Jj_i + hj (j = 1, 2, . . J) (71c) 

where k^, r^, and hj are defined in figure 8. . .. 

The difference equations which are to approximate equations (70c) to (70f) are 
obtained by averaging about the midpoint /X i,z. j\ 

I ""1 ^'2 ^' 2 j 

J jn,l ^ 

J = (72a) 

h, 


,n,i „n.i 


ui‘>‘ - uV’i 
i M 






(72b) 


g "’' - gM 
fL_^ = w.n,i 


^2 


(72c) 




'j 


H 


(72d) 


where, for example, 

„n,l lfcn,i + „n,l\ 
j-1 2lj J-lj 


The difference equations used to approximate equations (70a) and (70b) are rather 
lengthy. To illustrate the difference equations, an example equation similar to equa- 
tions (70a) and (70b) is chosen as follows: 


v’ + P jfv = x(u -^ + Pr7W — 'j 
\ 9X 9Z/ 

The difference equations for this equation are 


n-|_ 


V4 - Vi: i _ 

T ^ ^"5* I 2 


n-1 






ru, - u._i) 


(73) 


( 74 ) 


where, for example. 


V = i^p,i + 

J 4\3 J 3 3 / 
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Uji * 

"■” 1 = i (^>‘ * '■>-1 ^ *’"- i ) 

^“2 

The boundary conditions for equations (70) evaluated at x = and at z = zj are 


II 

0 

0 

II 

« 8 -‘ 

wS«i = 0 

II 

”S’‘ 


(75) 


j* /Vn-l,i-l „n-l,i-l „n-l,i-l „n-l,i-l tn-l,i-l\ /fn,i-l uP,i-l 

and ff?"^»\uP“^>Sv?"^»^,gP"^>^,wP“^»^,tP"^>^Vare assumed to be 
’ ] ’ j )• \ j ’ j ’ j ’ j ’ j J 

known for 0 ^ j ? then the difference equations (70a), (70b), (72), and (75) yield an 
Implicit nonlinear algebraic system of . 6J + 6 equations in as many unknowns 
^fj*,u?,v?,g?,wj^,tj’^ . This nonlinear system is solved by means of Newton’s method. The 

resulting linearized system is then solved very efficiently by using the block elimination 
method discussed by Isaacson and Keller (ref. 18). 

Numerical Formulation of the Energy Equation 
A new dependent variable G(x,z,? 7 ) is defined as 

G = E' ^ (76a) 

and equation (61) is written as 

uiE _G^ + p„(w|S -G|i)l (76b) 

8x ax ‘V 9z 3z/ 


(MiG)’ + M2G + Ms = xPio 


The difference equation for (76a) is written again by averaging about the midpoint 
^XjjjZj,?; and is similar to those given by equation (72). The difference equation 

for equation (76b) is written similar to equation (74). The boundary conditions for an 
adiabatic wall are 


gJ'^ = 0 


eJ»^ = 1 


(77) 
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The resulting algebraic system of 2J + 2 equations in as many unknowns 
which is linear, is directly solved by the block elimination method. 

RESULTS AND DISCUSSION 


One obvious difficulty in evaluating the accuracy of the three-dimensional turbulent 
boimdary -layer calculations on wings is the lack of complete, reliable data. Fortunately, 

. however , there are ajew good data available for flows. with simple geometries. Calcula- 

tions for these flows serve the useful purpose of evaluating the turbulence models used 
for the Reynolds stresses. References 3 and 4 present several comparisons of calcula- 
tions with experimental data. Although these comparisons are for flows over simple 
geometries and the calculations are for a coordinate system different than the one consi- 
dered here, the generally good agreement observed in those calculations gives some con- 
fidence in the accuracy of the turbulence model used in this study. Figure 9, taken from 
reference 19, shows the flow geometry and comparisons of calculated and experimental 
results for a planar turbulent boundary-layer flow approaching a three-dimensional 
obstacle. The results shown are for velocity profiles in a gradually steepening adverse 
pressure gradient flow off the plane of symmetry. The calculations were made for a 
Cartesian coordinate system which can be obtained from the present equations by setting 
hi = h 2 = 1, Ki = K 2 = Ki 2 = K 21 = 0, and 0 = Ti/l. 

Skin-friction coefficients are presented in figure 10 for the upper surface of a swept 
wing whose planform is given in reference 10. The calculations were made by obtaining 
the velocity components from the experimental pressure distribution by the procedure 
discussed earlier. To simulate the actual geometry, a reasonable thickness distribution 
was added to the planar wing considered in reference 10. As in reference 10, the calcu- 
lations were made for a unit Reynolds number of 4.92 x lO^ per meter and for a free- 
str earn Mach number M«> of 0.5.- In the figure z = 1.778 m represents an inboard 
station on the wing and z = 4.572 m represents a station in the middle of the outboard 
panel of the wing. The skin-friction coefficients are defined as surface shear-stress 
components normalized with free -stream dynamic pressure. Here, Cf^^ represents 
the shear-stress component in the x-coordinate direction and C{ ^ represents the shear- 
stress component normal to the x-coordinate in the tangent plane. In physical and trans- 
formed coordinates, they are defined by the following formulas: 


^,x “ 


Tjj + Tg cos 9 
ip U 2 

2 ^00 



UpWa M a 

+ gw ® 


Uoo 


(78) 


%z 


sin 9 




^^w /^e\^e^e 


g!^ sin 9 


(79) 
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Although the present results and those of reference 10 are qualitatively similar, 
there are several quantitative differences between the two predictions. One possible 
reason for the differences could be the starting procedure used to compute the initial 
conditions along the spanwise direction. Our calculations were made for a turbulent flow 
starting at approximately 3 -percent chord whereas the calculations of reference 10 were 
made with transition to turbulent flow occurring at 10-percent chord. ' Another possible 
reason for the differences could be the procedure used to get the velocity components 
from the experimental pressure distribution. 

The method described in the previous section resulted in the three computer pro- 
grams which are used separately from each other. One computer program deals with 
the calculation of the velocity components from the experimental pressure distribution by 
using the sweep theory. Obviously if the velocity components are known from the invis- 
cid flow theory, then this program is not heeded. The second computer program deals 
with the calculation of the nonorthogonal coordinate system and its geometric parameters, 
namely, the metric coefficients hj, h 2 , Kj, and K 2 appearing in the governing- 
boundary-layer equations. Through the use of this program, the coordinate system and 
its geometric parameters are calculated once and for, all for a given wing. The data is 
punched out on cards to be stored. If no changes are made in the airfoil cross sections, 
then this data can be used for any number of boundary -layer calculations without using 
the second computer program again. The third computer program deals with the solu- 
tion of the governing boundary-layer equations for a nonorthogonal system using the very 
efficient and accurate Cebeci -Keller Box method. This program assumes that initial 
conditions on two intersecting lines are given. In the present program, the two intersect- 
ing lines correspond to the wing -fuselage junction and to a line along the span a small 
fraction of the chord length away from the leading edge. This computer program solves 
the boundary-layer equations in a surprisingly small amount of time for a given external 
velocity distribution (either experimental or theoretical) and for a given wing coordinate 
system for both incompressible and compressible flows. The results in figure 10, for 
example, were obtained for a wing consisting of 29 z -stations and 19 x-stations with 
30 7/-points across the boundary layer. The. total central -processing -unit (CPU) time for 
all stations is approximately 30 sec oh an IBM 370/165 computer. 

FUTURE WORK 

The method described here has been tested for only one flow condition. It lacks 
certain important features and capabilities that may become very useful at different flow 
conditions, particularly for the third computer program which solves the boundary -layer 
equations. These features and capabilities conveniently can be divided into three separate 
tasks. 
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1. It is desirable and useful to include the canahiUty of starting the calculations at 
the stagnation line rather than some small distance aft of the stagnation line as in the 
present procedure. This task involves the solution of a special set of equations, called 
attachment -line equations. With this capability, the solution of one of the initial lines 
(sta^ation line) becomes exact but remains approximate (though a 'good approximation) 
on the other initial line (wing -fuselage juncture) as before. 

2. In the present method, the dimensionless cross -flow velocity is defined by 
g'-= w/we- However, this definition is not very convenient. In some problems where 
the outer- velocity component Wg changes sign," certain amibiguities arise. For exam- ' 
pie, if the cross flow at the outer edge of the boundary layer becomes slightly negative 
but remains positive in the rest of the boimdary layer, the value of g' will suddenly 
change sign from one station to the next. This introduces some discontinuity in the flow 
field since as Wg goes through zero, the value of g' becomes infinite at some net 
point between the two calculation stations. To avoid this problem, the transformation 
needs to be changed slightly and the cross -flow velocity w normalized by some refer- 
ence-velocity which does not change sign. 

3. A very important study that needs to be conducted involves the procedure with 
which the calculations are advanced in the spanwise direction. In the present program, 
a special solution at the root station is obtained prior to calculating the boundary layers 
on consecutive spanwise stations. At each spanwise station the solution starts with an 
initial profile and proceeds along the chord until We becomes negative. At that point, 
the program proceeds to the next spanwise station and initiates the calculation at the 
leading edge and so on. With this procedure the wing is covered from the root to the tip. 
It should be noted that region I is defined to be the region where Wg is positive. The 
calculations in region II (this corresponds to the region where We is negative) start 
from the -Wing tip.. The same approximate boundary-layer equations are solved as for 
the wing -root section to generate the initial conditions along the chord at the wing tip. 

The rest of the calculation procedure is identical to region I except that now marching is 
in the inboard direction as the boundary layer is calculated in consecutive spanwise sta- 
tions all the way to the wing root. 

This procedure of marching back and forth requires further study. If there is 
another region where the cross -flow velocity We changes sign, proper Ic^ic must be 
incorporated in the computer program. 

An alternative procedure to define separate regions can be utilized by the appear- 
ance of negative cross -flow velocity. In such cases, a procedure similar to the present 
marching procedure can be used. The proper marching procedure requires an extehsive 
and careful study since the locus of streamlines is unknown a priori on complex geome- 
tries. An efficient method can only be found by making the actual calculations, and chang- 
ing and testing the Ic^ic as required. ' 
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Figure 7.- Notation for conical section ABDC and 
the coordinate system in the developed plane. 
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Figure 8.- 'Net cube for the difference equations 
for three-dimensional flows. 













